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Abstract 

In this work we study the influence of populations mobility on the spread 
of a vector-borne disease. We focus on the chikungunya epidemic event that 
occurred in 2005-2006 on the Reunion Island, Indian Ocean, France, and 
validate our models with real epidemic data from the event. We propose 
a metapopulation model to represent both a high-resolution patch model 
of the island with realistic population densities and also mobility models for 
humans (based on real-motion data) and mosquitoes. In this metapopulation 
network, two models are coupled: one for the dynamics of the mosquito 
population and one for the transmission of the disease. A high-resolution 
numerical model is created out from real geographical, demographical and 
mobility data. The Island is modeled with an 18,000-nodes metapopulation 
network. Numerical results show the impact of the geographical environment 
and populations' mobility on the spread of the disease. The model is finally 
validated against real epidemic data from the Reunion event. 
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1. Introduction 



Many factors have influence the emergence and re-emergence of vector- 
borne diseases [D, 0|. Brutal changes in natural habitats so as massive or 
recurrent population migrations tend to speed up the spread of vector-borne 
diseases. The spatiotemporal evolution of such diseases is becoming a key 
issue for epidemiologists. To this purpose, models that consider the spatial 
distribution of a natural environment are of great interest. In this paper we 
are interested in the spatial spread of a vector-borne disease under the effect 
of human and vector mobility. Indeed, human migration is one of the factor 
that have influenced the re-emergence of several diseases 0, 0|. The mod- 
eling of geographical environments and populations mobilities are becoming 
mandatory in this context. There are several approaches to describe such 
spread. One of the typical approaches, introducing spatial spread variation 
in epidemic models, involves the use of partial differential equations 0,0, 0|. 
However, in the case of human mobility these approaches may not be appro- 
priate. The theory of "metapopulation", first introduced in 1969 [sj], in the 
field of ecology, allows such modeling. Several researches have been devoted 
to the study of disease spread in heterogenous environments 0, [kJ U|- In 
12l . Il3| . authors study the influence of human dispersal among n patches 
in the dynamics of disease spread. In |l4j, authors, to study the spread of 
seasonal influenza, focus on air displacements and model the environment 
with a network where nodes are airports and edges represent flights. In 15 



and |16|, authors are also interested in the case of influenza and mobility in 



the authors propose a virus spread 
the mobility of mosquitoes has been 

Other studies 

In 




terms of long journeys, 
model based on this theory. In 

modeled to study the spatiotemporal dynamics of malaria, 
focus on direct-transmission diseases like 



19, 17 



14 , influenza in the 



case of long trips (aircraft flights) is tackled. In 22] , the authors rely on the 
Ross-Macdonald model 
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to consider human mobility. 
Human mobility is usually considered in epidemic problems. Also, models 
that consider the dynamics of the vector population are numerous. However, 
to our knowledge, no model considers to couple population dynamics with 
populations mobility as we propose here. Moreover our approach promotes 
the consideration of vectors's mobility since the resolution of our model is 
greater than usual models. 

In this work we are proposing to couple two models (published in (25|): 
a mosquito dynamics model (growth and evolution of the population) with a 
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transmission virus model between two populations (humans and mosquitoes). 
The dynamics of the mosquito population is described by a stage structured 
model based on its biological life cycle. The different compartment states 
are egg, larva, pupa and adult. The disease transmission model is, for the 
human population, a SIR (Susceptible, Infected, Recovered) model. For adult 
mosquitoes, the transmission model is a SI (Susceptible, Infected) model. 

Those two models are formalized using the metapopulation theory. It 
considers a network where nodes represent real habitats of the environment. 
In each of these nodes, transmission and population dynamics models appear 
and are coupled with neighbor nodes. Links in the network represent both 
the local neighborhood of a node and farther nodes that code for the mobility 
of humans. 

We focus on a real case of chikungunya epidemic with the 2005-2006 event 
that occurred on the Reunion Island, a French island in the Indian Ocean. 
The island is modeled with a network. Since we want that network to reflect 
the local population's density, we consider the road network of the island as 
a proxy to the human density, considering that each crossroad is a node of 
the network. Then the local population on each node is adjusted according 
to real data given by the French Institute for Statistics (INSEE). Finally the 
all island is modeled with a 18,000 nodes network. 



In [26[, authors propose various distributions that match a data set of real 
human mobility patterns (cell phone probes). We rely on these distributions 
to create human mobility patterns for the population in our model. We 
assume that individuals only change disease status when they are on a node 



and not during displacements. See |27[ for a model with disease transmission 
during travel. 

Results show the decisive influence of mobilities over the spread of the 
disease. Not only the human mobility but also the vector local interaction 
that play an important role at the considered scale. Results are then com- 
pared to real data epidemic regarding the 2005-2006 event at the scale of the 
all Reunion Island and the model is validated. 

The remaining of this paper is organized as follows. Next section recalls 
original transmission and population dynamics models that this work is in- 
spired of. Section [3] introduces the original metapopulation model that is able 
to include the previous models in a network of patches, linked to represent 
populations mobility. Then, in Sect. HI according to the wish to validate the 
model against a real epidemic, a numerical implementation of the metapop- 
ulation model is constructed. This section details the construction of the 
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metapopulation network in terms of environment and populations mobility. 
Section presents various analysis performed on the model, showing the im- 
pact of mobilities. A validation is then proposed with a comparison with real 
epidemic data. Then a stochastic analysis of the model is proposed to study 
the robustness of the system. Finally Sect. |6] concludes that paper. 



2. Original Model 



We first recall the model proposed and studied in |25[ describing the 
vector dynamics. The formulated model is a stage structured model based 
on the biological mosquito life cycle. The vector population is subdivided into 
several classes: the aquatic stagese consisting in eggs E and larvae/pupae L, 
and then the adult stage A, representing adult females. We assume first 
that the number of eggs b, laid by females, is proportional to the number of 
females itself. Secondly, the number of eggs and larvae is regulated by the 
effect of a carrying capacity Ke and Kl, respectively. Then we formulate a 
transmission virus model where the adult female stage A is subdivided into 
two epidemiological states: susceptible S m and infectious I m . It assumes that 
there is no vertical transmission of the virus, so that births from susceptible 
and infectious mosquitos occur into the egg stage with the same rate. 

The human population consists in three epidemiological states: suscep- 
tible (or non- immune) Sh, infectious In, and removed (or immune) Rh- It 
is assumed that there is no vertical transmission of the disease, so that all 
births occur into the susceptible human class, at the rate bn > 0. Moreover 
we assume that the total human population Nh remains constant, thus birth 
and death rates are equal. An infected human is infectious during 1/jh days, 
called the viremic period, and then becomes resistant or immune. 

Forces of infections used in the model, which describe the rates of appari- 
tion of new infections, are standard and modeled by the mass-action principle 
normalized by the total population of humans, given by /3 m lH(t)S m (t) /N H 
and ^HlmifySnit) / Nh where (3 m and 0h are the transmission parameters. 

This hypothesis are summed up in Fig. [TJ 

Based on our model description (see Fig. [TJ and assumptions, we establish 
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Figure 1: Compartmental model for the dynamics of Aedes albopictus mosquitos and 
the virus transmission to human population. b(t) = bA(t) (1 — E(t)/ Ke), s(t) = 
sE(t) (1 — L(t)/Ki,) and the parameters of the model are given in table [1] 



dR H (t) 
dt 



the following equations. 

d -i {t) = bA(t)(l-^)-(s + d)E(t) 

§W = sE(t) !l-M\ isL + dL)m 
dA 

— it) = s L L(t) - d m A(t) 

^(t) = s L L(t) - d m S m (t) - p m ^S m (t) 
ctt l v n 

TT W = fi"rn~T7 S m (t) — d m I m (t) 

dt Nh 

^f(t) = -(5 H h^S H (t) + b H (S H (t) + I H (t) + R H {t)) - d H S H (t) 
dt Nh 

^ = Pu^Snit) -^u^-duluit) 



-yI H {t) - d H R H {t) 

(1) 
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Table 1: Parameters of the mosquito and transmission dynamics models. Most of the 



Parameters 


Description 


Value 


b 


oviposition rate 


6.0 


K E 


egg carrying capacity 


1000 


K l 


larva carrying capacity 


K E /2 


s 


transfer rate E — > L 


1/3 




transfer rate L —> A 


1/10 


d 


egg mortality rate 


1/3 




larva mortality rate 


1/3 




female adult mortality rate 


1/14 


b H 


human birth/death rate 


1/(78*365) 


Ph 


infection rate : vector to human 


See Sect. [53] 


P m 


infection rate : human to vector 


See Sect. EE 


1H 


human recovery rate 


1/7 



The study of the model is detailed in 25j . For this model (Op 
basic reproduction number is given by 



the second 



R 



d m {l + b H ) 



H 



1 



sK E s L K L 



d m (sK E + (s L + d L )K L ) 



where r = (b/(s + d))(s/(s L + d L ))(s L /d m ). 

The threshold r governs the dynamics of mosquitoes. In this article we 
assume that r > 1. This ensures the existence, persistence and global stabil- 
ity of a unique endemic equilibrium (E*, L*, A*), which corresponds to the 
biological and interesting case. The second reproduction number Rq governs 
the dynamic of the transmission model. 



3. Metapopulation Model 



In [18|, [19|, [17| J. Arino et al. formulate a general system of differential 



equations allowing to describe human mobility. In this model, they identify 
each population by its origin and its present location. We rely on this model 
and extend it with the definition of a neighborhood in which humans and 
mosquitoes interact. 

3.1. Human mobility 



Based on the model given in [18|, [19|, [17(, we proposed an extension of our 



previous model [251 ]. to describe the spread of a chikungunya under human 
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and vector mobility on a large scale network. 

Assume that the number of nodes in the network is n. A human pop- 
ulation is identified due to two characteristics: the node from which he is 
originated, i. e. its residence and the node where he is at time t. We assume 
that the total human population is constant, i.e. births and deaths occur 
with the same rate bu = du- Moreover, we suppose that birth occur in the 
resident node while deaths take place in any node the human is present. 

Let us denote by Smj(t), Iuijif) and Rmjif) the susceptible, infected and 
removed human populations originated from node i and present on node j 
at time t. We denote by S mi and lu% the susceptible and infected mosquito 
present in node i. In a first assumption mosquitoes mobility is neglected, 
which looks like realistic compared to the distance of humans displacements. 
So, resident mosquitoes from node i are also present on this node. 

The total number of susceptible, infected and removed human residents 
of node i, is given respectively by S H l = J2]=i S Hij, IhI = Y^j=i I Hij and 

D r D 

K Hi — 2-^=1 K Hij- 

The total number of susceptible, infected and removed humans present on 
node i, is given respectively by SW? = YTj=i S Hji, Im = YTj=\ I Hji and 

The number of human residents on node i is then equal to Nh[ = ShI+Ih^-'t 
RjiI and the human population present on node i is N H ^ = S #f + 1 #f ' + RhI '. 



As in [32( and [17], we define the travel rate from node i to node j by 
girriji, where gt > corresponds to the per capita rate at which residents of 
node i leave this node and a fraction rriji > of them go to node j, with 
ma = 0. Residents of node i, present on node j, then return to node i with 
a per capita rate > 0, with ru = 0. Displacements between two nodes are 
represented in Fig. 

The dynamics of human populations (susceptible, infected and removed) 
resident and present on node i is given by: 

77 ~ = d H (N H \ — Shu) ~ 9%Sm% + r ikSmk — ^Hi^rr^SHu 

tt = ~d H lHu — gdua + r ik I Hik + (3 Hi i Shu ~ IrIru 
dt <H N H P i 

— 77^" = —dHRuu + JhIhu — 9iRhu + r ik ^ H ^ 

k=l 
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Node i 




Figure 2: Human population mobility between nodes i and j. 

The dynamics of human populations (susceptible, infected and removed) 
resident on node i and present on node j is given by: 



dS H 



dt 
dlmj 

Hij 



dt 
dR 



—dnSnij + giTTijiSHu — t ijSnij — fiuj N p Snij 

j 

Itti ' 



dt 



—diiR-Ha + //,"',, /«'//,, — r ijRmj + IhIh 



The dynamic of mosquito populations (eggs, larvae and pupae, adult 
females) given in 25| in each node is: 



dSmi S mi 

at l\ Hi 

dl-mi 



j P j j 

mi N p u "m 1 mi 



Hi 



dt 

b(S mi (t) + I mi (t)) (l ~ ^ ) ~{s + d)E t {t) 



dt \ K E 



(2a) 
(2b) 
(2c) 
(2d) 



where the dynamic of immature stages E and L is described by corresponding 
equations in system (Tj[|). 

Remark 1. i. Note that the infection transmitted in node i between sus- 
ceptible mosquitoes and infected humans depends now on the human pop- 

n S S 
ulation present on this node Nh\: / Pmi J^p ^Hji = firm J^p IiU- 

~[ N Hi N Hi 
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ii. The distance between nodes is not explicitly taken into account, never- 
theless it is implicitly included in the coefficient rriji and . 

hi. In the following we focus on daily displacements, humans who leave their 
resident node, obviously return to their resident node on a daily basis. 
Displacement matrices M T = [giiriji] and R = [r^] have the the same 
zero /nonzero pattern. 

The human mobility model and the virus transmission dynamics is given 
by system fl3]). 

j, U =dH(N H r i—SHn) —gjSHii + 'y^ J f'ikSHik — ^Hi-^-^SHii (3a) 



-T~ —9i m jiSHii — dHSHij — TijSuij — 0Hj AT ^p Suij (3b) 

at iv h 



"d 



dI H - n Im- 

—— =—dHlHu—9ilHii+/ r ik I Hik +(3 Hi p Shu—JhIhu (3c) 
dt N H l 



^ —9i^fjHii — dHlHij—rijlHij+PHj^-^pSHij—lHlHij (3d) 



I'J 'IJ^-Ulj I HHj A j j^tl % ] 

"j 



dR H - n 

— ~t^~ = IhIhu — duR-Ha — giRun + ^ r ik R Hik (3e) 

fc=i 

dRuij 

— ^ — = gi^njiRHii + lHlHij — dHRHij—rijRHij (3f) 

3=1 



mi j * , ^ \ ^ n °mi 1 Hji /„ x 

—SlLh —a m Dmj — / J Pmj j^^p l<JgJ 



/ j Hm% *j p u 'm 1 mi \ OLL ) 



dt /-^ ' 1 Nj 

3=1 

i r * A 1\ ^Li , 1 . (s L + d L )K Li . . 

where L* = 1 and 7^ = 1 H is given by the 

V r JlLi sK Ei 
endemic equilibrium of the vector population dynamic model. In this model, 
each node is described by (3n + 2)n equations. 

3.1.1. Equilibrium of the model 

Proposition 3.1. The nonnegative orthant ]R + ( 3n+2 ) n is positively invariant 
under the flow of (j3]) and, for all t > 0, Shu > and Smj > provided that 
girriji > 0. Moreover, solutions of (|3]) are bounded. 
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Proof. We easily see that solutions of system flS]) remain nonnegative. Indeed, 
it is sufficient to show that for all nonnegative initial condition, the vector- 
field points out to the interior of the positive orthant. Assume now that 
Shu = at t = 0, then dS H u/dt = g^A^ + XILi r ikSuik > 0, thus S H u > 
for t > 0. Equally, if Shh = at time t — 0, then dSmj/dt = gimjiSun > 0. 

Finally, the boundedness follows from the positive invariance of R + ( 3n+2 ) n 
and the constant population property Indeed dN^l/dt = 0, which means 
that, for any node i, the resident population is constant, thus by extension, 
the whole population is constant. Moreover, solutions of ()3]) are bounded, 
since (M+)( 3n+2 ) n is invariant, human population is constant and vector pop- 



ulation is bounded 25 . □ 



Definition 3.2. tla l 1. The system is at an equilibrium if the time deriva- 
tives in ([3]) are zero. 

2. A node i is at the disease free equilibrium (DFE) if Iuji — 0, I m i — for 
all j = 1, n. 

3. The n-nodes model given by ([3]) is at the DFE if each node is at the DFE, 
i.e., Inji = and I mi = 0, for all i,j = 1, ...,n. 

Proposition 3.3. System ([3]) always has the following disease free equilib- 
rium: 





-(• 


j * 


= o, 




= o, 




SL 
dm 




NhI, SH*j — gi- — — — ShI 



rriji 
d H + r { 



L*, lm* = 0, for alii, j = l,...,n, j. 
Proof. It is sufficient to remark that 

-i v^'« m ki r ik v-^™ \-^n m ki r ik , sr^n mn r— i 



Theorem 3.4. }laj . Assume that system (jHJ) is at an equilibrium and a 
node i is at the DFE. Then all nodes that can be accessed from node i and all 
nodes that have an access to node i are at the DFE. Moreover, if the outgoing 
matrix M T is irreducible, then all nodes are at the DFE. 

Proof. See appendix. □ 



10 



Definition 3.5. lla l. The disease is endemic within the population if the 
number of infective individuals is positive in this population. 

The disease is endemic on node i if there is a population on node i in which 
the disease is endemic, i.e., there exists k £ {1, ...,n} such as Ihh > 0. 



Theorem 3.6. llal . Assume that system §3§ is at an equilibrium and the 
disease is endemic on node i. Then the disease is endemic on all nodes that 
can be reached from node i. In particular, if the matrix M T is irreducible, 
then the disease is endemic in all nodes. 

Proof. See appendix. □ 

3.2. The vector mobility model 

Contrary to human displacements, it is unrealistic to identify mosquitoes 
by their origin and destination. Nevertheless, we know that Aedes albopictus 
mosquitoes have a limited flight range. Most mosquitoes disperse less than 



two hundred meters away from their original breeding place |33|, |34| . Indeed, 
mosquitoes present in a node i may have an activity in neighboring nodes j, 
depending on their proximity (defined later). Metapopulation models of the 



spread of vector-borne diseases (see for instance [35|, |36|) that include human 
long-distance displacements, do not take into account mosquito mobility. In 
our case of daily movements and very precise resolution, we cannot neglect 
the influence of mosquitoes displacements. We propose to model this activity 
using the biological radius of interaction around their breeding sites. Let us 
denote by the distance between nodes % and j and by d max the maximum 
interaction radius of mosquitoes (approximately 200 m). In particular, we 
have da = for all % = 1, ...,n. Now we assume that mosquitoes originated 
from node % interact with population of node j, according to a function of 
the distance linearly decreasing. This function is given by 

{ dmax ~ dij 
dmax (4) 
else 

Then, model (El) becomes 
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which describes both the vector population dynamic and the virus transmis- 
sion, through the human and vector mobilities in the network. Note that 
aquatic stages of eggs and larvae remain described by equations (l2c|) and 
fl2d|) and complete the model. 

4. Application 

Willing to validate this approach, this section proposes a comparison of 
this model with a real-life example of chikungunya epidemic, namely, the 
event that occurred in 2005-2006 on the Reunion Island, Indian Ocean. This 
validation process makes a strong usage of real and realistic data, to reflect 
as much as possible the original environment. This last one is modeled based 
on real geographical information. Populations mobilityis modeled following 
existing eal data analysis. Finnally, epidemiological results are compared to 
real data from the 2005-2006 event. 
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4-1. Distribution of the Human Population 

The Reunion Island is a mountainous region and the local population is 
not uniformly spread over the land. To reflect this particular density in the 
metapopulation model, we rely on two realistic sources of information. 

The INSEE gives access to an estimated density of the population, based 



on a mesh zoning of the space |37|. This zoning consists in squares of one 
kilometer long. In each square, the number of people living in is estimated 
based on both tax cards and cadastral data. Figure [3] shows these 1 km 2 
squares with their associated population. 



P V-- 



1 to 19 ^ y 

20 to 199 i 
200 to 499 
500 to 999 
1000 to 1499 
1 500 to 4999 
5000 to 14999 
>15000 



■ 
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Figure 3: Estimated population distribution of the Reunion in 2007, according to the 
French Institution for Statistics, on a 1 km 2 granularity. 

There is a good confidence in the quality of this data, however, the limited 
granularity of one square kilometer is not precise enough for our purpose. 
Indeed this model envisions a lower scale of interaction between humans and 
mosquitoes, since the latter only have a couple hundreds of meters interaction 
range. 

Secondly, the road network is used as a proxy to the human density. 
Indeed, the road network is denser where the population is itself dense, and 
lighter where no one lives. Data extracted from OpenStreetMap (OSM) (38| 
provide road information and especially road intersections at a high resolution 
(scale of one meter). We propose to use those intersections as the nodes of 
the metapopulation model. 
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Then the human population is distributed on these nodes. Each of them 
belongs to one of the 1 fcm 2 -areas of the INSEE. So, on each of these squares, 
the local population can be evenly distributed among its own nodes. If a 
square has no node associated to it (it happens on low density areas) then a 
single node is created with the corresponding population. 

Finally, this distribution is at least as relevant as the real data given by 
the INSEE at the scale of 1 km 2 , but the resolution goes below the 1 km 
limit, thanks to the crossroad/node analogy. Figure H] illustrates the main 
steps of the construction of the model. 




(a) 



(b) 



Figure 4: Schematic construction of the nodes of the metapopulation network. 4(a)| 
network nodes are the roads intersections, knowing that this network may not be totally 
accurate (red roads missing in the map). |4(b)| distribution of the population on nodes 
according to the INSEE data. Empty cells (with no known road network) are given one 
node (as in Fig. |4(b)| second row, rightmost cell). 



4-1.1. Distribution of the Mosquito Population 

No information about the density of mosquito population is available. 
Moreover, we are not interested in the whole population but rather in the 
part of it that interacts with humans. Our hypothesis is that the density of 
mosquitoes is homogeneous. 

To retrieve this density from the human-centered network, each node is 
given a number of mosquitoes (through the carrying capacities) proportional 
to the geographical surface associated with that node, thanks to a simple 
Voronoi tessellation. 

Since we only consider mosquitoes that interact with humans, the surface 
is upper-bounded with the disk of radius d max that expresses the maximum 
interaction distance for a mosquito along its life. 
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Let Si be the surface of node i. Let S max = ~K(ft max be the surface associ- 
ated with the maximal interaction radius d max . Carrying capacities for each 
node are then set as follows: = Ke <f(Si) and K^ = Kl <f(Si), with 
Ke and Kl, constants for the carrying capacities of eggs, respectively larvae 
and ip(Si) is equal to Si/S max and is upper bounded by 1. 

4-1.2. Human Mobility 

In the model, humans trips are given by outgoing matrix M T and incom- 
ing matrix R and represent for any pair i and j, the probability for humans 
living on i to go to j. Non-null values in the matrix define edges between 
nodes in the metapopulation network. 

To reach realistic motion we need to rely on realistic data or model for 
human mobility. Unfortunately, we could not get such information for the 
Island, so we rely on some more general analysis work from M. C. Gonzalez, 



C. A. Hidalgo and A. Barabasi 26] who analyze mobile phone communication 
logs of 100,000 users. These logs register the geographical position of mobile 
phones when calls happen or when texts are emitted/received. Mobile phones 
are a good proxy to the human mobility because users always carry their 
phone. 

The authors were able to produce general laws on observed mobility pat- 
terns. These results give general formulas that describe the mobility ob- 
served. We use these formulas to generate new human mobility patterns. 
We consider the three following measures. 

Trips length. The distribution of the length of a human trip Ar in kilometer 
according to the following power law: 

P(Ar) = (Ar + Ar )" /3 exp(-Ar/K) (6) 

where Ar = 1, 5km is the cutoff value of the law and k = 80km. 

Presence probability. Given N possible destinations for each individual, 
their presence probability in each destination is approximated with the fol- 
lowing Zipf law 

f( k]N ) = —^ 

Eti(Vn) 

where k is the rank of each destination when decreasingly sorted. 
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Return probability. The authors observe that there is a peak of probability 
of returning to the same place every 24h hours. In other words, the human 
mobility mainly follows a daily pattern such as home/work trips. 

From these measures we generate artificial per-individual trips on the 
island that respect the observed properties in the original data set. 

Remark 2. The model defines matrices M T and R as tables indicating re- 
spectively departures and returns. The mobility model we could create out 



from \2b\ l however gives the probability of presence for a human, given all of 
its possible destinations. Our hypothesis here is that there is a link between 
departures / returns probabilities and presence probabilities that we do not 
investigate in this paper. 

4.1.3. Mosquito Mobility 

As stated in the model description, mosquitoes are not identified with 
origins and destinations. Since we are aware of their limited flight range, their 
mobility is thus defined by a geographical disk area of interaction centered 
at their origin node. Function fll]) defines their interaction with humans. 

This new mobility pattern is constructed, in the metapopulation network, 
with edges linking nodes that are below the interaction range d max . 



5. Results 

In this section we show and discuss results of the simulation of our 
metapopulation model applied to the scenario of the Reunion Island de- 
scribed in the previous section. Our results are then compared with real 
data from the 2005-2006 epidemic event. 

5.1. Analysis of the Metapopulation Network 

The metapopulation network can be represented by two graphs sharing 
the same set of nodes, or in other words, it is a graph with 2 subsets of 
edges. One subset is for human mobility and one is for mosquito mobility. 
The purpose of Table |2] is to give some metrics for those two graphs to give a 
general overview of theirs dimension. Indeed due to its size no visualization 
of the all network is proposed since it does not give any useful information. 
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Tabic 2: Metrics for the two graphs of the mctapopulation network. 



Mosquitoes 



Humans 



number of nodes 
number of links 
average degree 
connectivity 
diameter 



no (1729 connected components) 
71 (for the biggest component) 



17988 
151772 
« 17 



17988 
744313 
« 83 
yes 
15 



5.2. Spread of the Disease in the Network 

In this scenario we observe the spread of the disease with one infected 
individual put on one node of the network. Provided infection parameters 
(Ph and j3 m ) are set high enough, the insertion of this individual in a system 
at disease-free equilibrium may start an epidemic event. Three nodes of the 
network are monitored: the first where the individual was inserted, a one- 
hop neighbor, and a farther one located 6km from the insertion node. Four 
metrics are observed: susceptible humans (S#), susceptible mosquitoes (S m ), 
infected humans (/#), and infected mosquitoes (I m ). 

Figure [5] shows, as expected, a shift in time of the evolution is observed 
from the closest nodes to the farther one. However, that shift (and thus 
the spread of the disease) is not proportional to the geographical distance 
between nodes but rather to a graph distance. The third observed node is 
60 hops far from the insertion node in the mosquito mobility graph and only 
4 hops away in the human mobility graph. Human mobility, thus, greatly 
speeds up the spread. 

5.3. Consequences of the Mosquito Mobility 

As stated in the introduction, most of the metapopulation models on 
vector-borne viruses focus on long distance and long duration journeys for 
the human population. In those models the vector mobility can easily be 
ignored. Here we consider daily-based mobilitty with short journeys. We 
believe this shorter scale in time and space brings the new constraint that 
mosquito mobility cannot be neglected anymore. 

The effect of local mosquito interaction is shown with a scenario where 
again, a disease free population is inserted an infected human. Two ex- 
periments are carried out, one with mosquito mobility and one without it. 
Human mobility is kept in both scenarios. 

Figure [6] shows the number of instantaneous infected humans (In)- The 
infection starts rapidly and is concentrated with mosquito mobility enabled. 
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Susceptible Humans(S H ) Susceptible Mosquitoes (S M ) 




25 days 50 days 25 days 50 days 



Figure 5: Spread of the disease in the network for three observed nodes. The first gets the 
infected human. The second is an immediate neighbor (125m away from the first). The 
third, 6km away from the first node, is 60 hops away in the mosquito mobility graph and 4 
hops away in the human mobility graph. Quantities of population (y-axis) are normalized. 
A shift (or delay) effect in the spread is observed. 

Oppositely, without this mobility, the spread takes more time to start and 
the peak is less important. It is clear that mosquito mobility at this scale 
has to be considered. 

5.4- Consequences of the Human Mobility 

Human mobility has an obvious effect on the mobility. The purpose here 
is to analyze realistic and slight modifications of this mobility like quarantine 
measures that may be taken in case of epidemics. 

This scenario proposes to stop the human mobility on infected parcels. So 
quarantine measures will be localized only on nodes where a given threshold 
of infection is reached. In quarantined nodes inhabitants are not allowed to 
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Figure 6: Effect of mosquito mobility of the spread of an epidemic. Evolution of the 
number of infected humans (Ih) with and without mosquito mobility. In order to obtain 
comparable results, infection rates as set to a relatively high level: (3h = 0.2 /3 m = 0.15. 

move out and no foreigners are allowed in. Since mosquitoes are not stopped 
by quarantine measures, they continue to interact within all nodes. 

Figure [7] shows global instantaneous and cumulated infections values for 
humans at the scale of the all population. In this scenario, the disease would, 
without any control, reach 35% of the population, just like the chikungunya 
event of 2005-2006. 

Results show the possibility to rapidly reduce the instantaneous infection 
rate. For instance, the Ih peak is almost divided by two when the quarantine 
threshold is set to 10%. However looking at cumulated infection cases (sero- 
prevalence), such a threshold does not significantly reduces the total amount 
of infected. To expect a real effect on the total number of infections one have 
to set the threshold below 1%. Here, each node on average has less than 50 
humans, so any threshold below 2% is technically impossible to achieve. 

5.5. Parameters Analysis 

The identification of predominant parameters of the system is needed 
to understand its behavior. We propose here a numerical analysis of the 
parameters of the system. An in depth statistical sensitivity analysis or the 
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Infected Humans (l H ) 




Start of infection 6 months 12 months 18 months 

Figure 7: Consequences of the human mobility. The instantaneous number of infected 
humans and the seroprevalence (or number of people ever infected) are observed during 
an epidemic. Infection values depend on the quarantine threshold that reduce, node per 
node, the human mobility, over the time. A threshold of 10% indicates that a node with an 
infection rate of at least 10% will be quarantined (incoming and outgoing human mobility 
are blocked). 

analytical estimation of parameters would be of great interest too, but are 
out of the scope of this paper. 

Infection rates parameters, /3h (mosquitoes to humans) and /3 m (humans 
to mosquitoes) have a strong impact on the results. Depending on the values 
selected for fin and f3 m , the infection goes from some few isolated cases to 
an epidemic that contaminates the whole population. 

Fig. El shows a scatter plot of values taken by (3 H and (3 m . For each 
couple ((3h ,(3 m ) the seroprevalence (after 400 days of epidemic) is observed. 
Values in percentage represent the ratio of the total human population ever 
infected depending on the two parameters. Contour lines help finding, given 
a seroprevalence percentage, values for /?# and f3 m . These vales are used in 
the validation. 
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Figure 8: Effect of the (3h and /3 m parameters on the seroprevalence. Blue crosses are 
values obtained by simulation. Contour lines are obtained with a bilinear interpolation of 
these data points. 



5.6. Validation against real data 

We propose a validation with a comparison against real data from the 
2005-2006 chikungunya epidemic. The two first cases of chikungunya where 
reported at the beginning of March 2005. The disease propagated during the 
following weeks and reached a peak at mid of May. The number of new cases 
then started to decrease until authorities thought the event was over at the 
end of the year 2005 with a seroprevalence (total number of cases) of 6,000 
people. But in the second half of December 2005 the spread started over with 
a strength that was not comparable to the previous peak. This second event 
reached a peak in February 2006 with more than 47,000 cases in a week. This 
sudden reactivation of the spread was later explained by a genetic mutation 



of the virus that lead to a new and more infective strain [39[. After the peak, 
the number of cases slowly started to decrease. The epidemic was declared 
over by April 2006. In the end, the InVS (French Institute for Health Care) 
counted 265,733 cases of chikungunya from March 2005 to April 2006. This 
represents more than 35% of the total population of the Island. 

The following results are compared to real data indicating, week per week, 
new cases of the disease. This information was kindly provided by the InVS. 

Back to our model, the new virus strain with its stronger infection rate 
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is modeled thanks to infection rates /3h and /3 m . We try to reproduce the all 
event by starting the simulation of the epidemic with one set of parameters, 
and then, by changing only one time the values of Ph and /3 m at the moment 
it appears in the real event. The two sets of parameters where selected 
experimentally with the help of the previous study (see Fig. [8]) giving possible 
values for (3 H and (3 m for a given seropre valence. 

Figure M compares the evolution of the real seroprevalence of the chikun- 
gunya event, with simulation results obtained with our model. Although the 
two curves do not fit perfectly, it allows a visual overall validation of the 
model. 



40% -] Real Data (Seroprevalence) 

■ - Simulation 



20% 



0% 




Figure 9: Comparison of the model with real data. The seroprevalence is considered for 
the duration of the Reunion Island epidemic. Values of @h and /3 m are 0.0118 and 0.0101, 
before the mutation event and 0.0245 and 0.0161, after. Those two couples (/3y,/3 m ) arc 
chosen thanks to the study in section 15.51 Values for the other parameters of the system 
are given in Table [TJ 

These results are a first step in the validation of this approach. In these 
experiments, only infection rates where investigated. Of course, other param- 
eters need to be considered. Especially those leading the human mobility. 
Moreover the metapopulation approach needs to be validated at a lower scale 
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than the scale of the all Island. These last observations form the perspectives 
of this work. 



5. 7. Stochastic Analysis of the Model 

The model can be affected by stochastic changes. Indeed, the human 
mobility given by matrices M T and R are defined according to the mobility 



probability law given in [26 



We propose to study the robustness of the system by running simulations 
of the model with various instances of the human mobility matrices and thus 
observe the variations in the obtain results. The population itself and the 
number of displacements do not change, only the destinations are different. 

Figure [TU] illustrates the evolution of infected humans on various (ran- 
domly chosen) nodes in the network, within the very same scenario that the 
mutation experiment described in section 15.61 We assume that observations 
made on In also apply on the other outputs of the model. 30 different human 
mobility matrices are used. The figure shows that the height of the variation, 
although different on each node, remains within coherent boundaries. Indeed 
when considering, for each node, the date when the standard deviation is the 
highest, compared to the node's population, it remains below 2% of that 
population, as illustrates Table |3j 



Tabic 3: Highest values of Standard Deviation (SD) for Ijj on a random selection of nodes. 





Node 1 


Node 2 


Node 3 


Node 4 


Node 5 


Population 


24 


16 


932 


21 


78 


SD 


0.27 


0.16 


0.76 


0.17 


0.13 


date of SD (in 2006) 


01-24 


01-16 


02-15 


01-12 


02-21 


% of population 


1.11 % 


1.05 % 


0.08 % 


0.84 % 


0.16 % 



When considering simulation results at the scale of the whole network, 
the variation is even lower with a maximum standard deviation (on the 10 4h 
of January 2016) of 544 infected humans which represents less that 0.07 % 
of the population of the Island. 

These results tend to confirm the robustness of the system against stochas- 
tic variations. 



6. Conclusion 

In this work the main concern was to try to validate dynamical systems 
for the spread of a vectorial disease against a real epidemic scenario. We 
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# / of <f* # ^ of <f° ^ 

Date 

Figure 10: Statistics on Ih for a random selection of nodes. The blue ribbon shows 
minimum and maximum values over the 30 experiments. Darker blue curves show 1 st and 
3 th quartiles while the red curve shows the median. 

focused on the Reunion Island's epidemic that occurred in 2005-2006. 

Many issues appear when trying to bridge the gap between global models 
and real life problems. Among them we chose to focus on the modeling of 
populations' mobilities. This logically implied a realistic modeling of envi- 
ronment, bringing us from temporal modeling to spatio-temporal modeling. 

Designed to consider spatial interactions, the theory of metapopulations 
allows the special geographical environment of the Island to be included into 
the original model. This metapopulation network was created based on real 
geographical and demographical data. It could then handle local interactions 
between nodes and thus allow populations mobilities to become part of the 
model too. 

We proposed two mobility models. The first one, for humans, is based on 
the analysis of real human mobility data sets (mobile phones probes). The 
second model, for mosquitoes, is based on local interactions between nodes. 
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After various studies of the consequences of mobilities on the spread of the 
disease, and after an analysis of the parameters of the system, a validation 
of the model against real seroprevalence data from the Reunion epidemic 
was proposed. Results validate the approach and clearly identify the human 
mobility as a key parameter in the spread of such an epidemic. 

As stated above, the choice in this work was to focus on mobilities. How- 
ever, other recognized issues play a key role in the all process. For instance, 
the effect of rain falls and weather seasons are known to directly influence 
mosquitoes evolution stages (especially aquatic phases). Studying the effect 
of seasons on the epidemic is definitely a perspective of this work. 

Finally, the integration of the metapopulation model and mobilities con- 
tributed to increase the complexity of the model. This last model which has 
only been explored on a numerical bases would need an analytical study. 
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Appendix A. 

Let V^ = {k^i : g(m ki > 0} and V-n = {k ^ i / g k m ik > 0} 

sets of nodes that can be accessed directly from city i and nodes that have a 

direct travel access to city i, respectively. Let 

Ai^ = {k ^ i/3(ki, ...,k g ) distincts, m kli m k2kl . ..m kkq >0 and gig kl . . .gk q >0} 
A^i = {k ^ i/5(ki, ...,k q ) distincts, m klk m k2kl . ..m ikq >0 and g k g kl . . .gk q >0} 

sets of nodes that can be accessed from city i and nodes that have an access 
to city i, respectively. 

Proof of theorem \3-4\ Assume that node 1 is at the DFE (without loss of 
generality), i.e. Iuki — 0, Wk = l, ...,n and I ml = 0. From eq. (I3cl) . we have 
dlji n /dt = Y^k=i r iklHik- But node i — 1 is at the DFE, i.e. dlnn/dt = 0, 
and since r Xv > 0, Wv G V]_>., then, I Hlv = 0, Wv G V\^,. Consider now 
equation (l3d|) with i = 1 and let v G Vi^, then 

dlHlv _ n I my Q 

dt N H P V 

As the system ([3]) is at an equilibrium point, thus dlniv/dt = 0. However 
(3hi > and Shi v > from proposition 13.11 then Im v = for all v G Vi_>. 
Finally, we have to show that for all v G and for all k = 1, n, iHkv = 
i.e., all humans resident of node k and present at node v are not infected. 
Consider equation (|3hl) for a node v G Vx_>.,then, 



dlmv n S mv I H j v it a ^mv j n ^mv j p p, 

"77 / j Pmv iy p Um-Lmv Pm v y p / j -^Hjv Pmv at p*H v 



Since (3 mv > 0, and from proposition ^. 11 S mv > then lu v v = 0, z. e. X)*=i -^a™ 

0. It follows that for all k = l,...,n, Ihhv = 0> since > for all 

1, j = 1, ...,n. Thus, all adjacent nodes to node % — 1 are the DFE equilib- 
rium. Moreover, by induction, we obtain that all nodes v G A±-> are at the 
DFE. 

Assume now that M T is irreducible. From eq. (13dl) with j = 1, we have 
dlnn/dt = giTriiilHa- Since system is at the equilibrium dlnn/dt = 0, beside 
girriu > for i G V^i, then In a = for all i G V^i. Let v G V-»i, from 
equation fl3cl) . we have 



rvklHvk + $ H i~T7~~p~ S f{ vv = 0, 
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n j. 

then \ r v klu v k — and (5 Hi mv p Sh V v = and Sh vv > 0, from proposition 

13.1^ thus L mv = for all v G V^.±. Finally, consider equation fl3h|l for v G V-^i, 
we have 



dl <? 

Ul± mv n u mv j p p, 

~dT ~ Pmv l^ v Hv ~ ' 

n 

then In* = 0, i.e. Ihjuv = VA; = 1, n. Thus all nodes are at the DFE. 

k=l 

Proof of theorem \3.6\ . 

Proof. Assume that the disease is endemic in node i = 1, i.e. there exists 
q G 1, n such that lu q \ > 0. We have to show that, in this case, Inn > 0. 
If q = 1 then we can proceed. Assume that q ^ 1. Assume by contradiction 
that Inn — 0- Since the system is at the equilibrium, from fl3c|) . we have: 

U — — -j— — + / r lk^Hik + PHi-rj—pdHii- 
dt ^-^ N H K 

k=i Hl 

Beside fim > and Shu > then I ml = 0. Thus, from equation (l3h~|) . we 
have 

_ (Urnl _ Y^/o g 

dt 2^p^ Nh v^i- 

Beside (3 mi > and S mi > then I H j\ = for all j = 1, ...n, which is a 
contradiction. We have Inn > 0, if the disease is endemic in node 1, which 
we now assume. 

Consider equation (j3dj) with i — 1 and j ^ i. Assume = 0. Since 
the system is at equilibrium, we have 



dl *n - nrri T IlTlj 
^ t — 9\ m j\ L HU + PHjJj—p 



If j G then g\vn^\ > thus Ihu = 0, which is a contradiction. Finally 
Imj > for all j G Vi_>, which means that the disease is endemic. Particu- 
larly, we deduce that Injj > from the first part of the proof. By continuing, 
we can show that the disease is endemic in all nodes j G A\^, that is to say, 
nodes reachable from node 1. □ 
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